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SECTr.ON  I 


TN7’ROrifCTIOK 


lYedictions  of  gau  turbine  comV'.uGtor  perl'crmnee  and  poilutar.t  enisoion 
characteristics  require  modeling  proced’ii’es  porscssing  a high  degi'ce  of  sophis- 
tication. fast  attempts  at  modeling  combustion  systems  have  been  largely 
frustrated  by  the  complexity  of  the  toupled  hydro'pynaLd.c  and  chemical  processes. 
The  difficalty  can  be  largely  attributed  to  the  lack  of  understanding  of  the 
flow  processes  which,  through  the  excliange  of  heat,  rmss  and  momentum,  con 
directly  relate  to  pollutant  forma  "-ion  and  combus‘:ior.  efficiency.  For  example, 
swirling  flow  has  been  shown  to  have  /in  j_mpoi‘tant  influence  on  the  stability 
and  combustion  intensity  of  flames  (def.  l)  as  well  as  on  residence -time  dis- 
tributions (Ref.  2)  in  combustors  vmich,  in  turn,  can  be  related  to  combustor 
performance  and  efficiency  as  well  as  to  emission  characteristics  (Refs.  3 and 
U).  Techniques  employed  in  modeling  combustor  flow  processes  have  generally 
been  highly  simplified,  particularly  flow  modeling  techniques  where  stirred 
reactor  concepts  and  one-dimensional  assumptions  are  employed  (Refs.  5 through 


1.  Beer.  J,  M.  and  N.  A.  Chlgler:  liability  and  Combustion  Intensity  of 

Pulverized  Goal  Flames  - Effect  of  Swirl  and  Impingement.  Journal  of  the 
Institute  of  Fuel,  December  1969- 

2.  Beer,  J.  M.  and  W.  Leucker:  Turbulent  Flames  in  Rotating  Flow  Systems. 

F^per  No.  Inst.  F-riAFTC-7,  North  American  Fuel  Technology  Conference, 
Ottawa,  Canada,  (1970). 

3.  Beer,  J.  K.  and  J,  B.  Lee:  'OiC  Effects  of  Residence  Time  Distribution  on 

the  Performance  and  Efficiency  of  Comoi^stors.  'The  Combustion  Institute, 

pp.  1187-1202,  (1965). 

U,  Marteney,  P,  J. : Analytical  Study  of  the  mnetics  of  Formation  of 

Nitrogen  Oxide  in  ifydrocarbon  - Mr  Combustion.  Combustion  Science  and 
Technology,  Vol.  1,  pp.  37-45,  (1970), 

5.  Fletcher,  R.  S,  and  J.  B.  HeywooU:  A Model  for  Nitric  Oxide  Emission  From 

Aircraft  Gas  Turbine  Engines.  A]AA  Paper  71-123,  (1971). 

6.  Hammond,  D.  C.,  Jr.  and  A.  M,  Mellor;  /'jialytlcal  Predictions  of  Emissions 
from  and  Within  an  Allison  J-33  Combustor,  Combustion  Science  and  Tech- 
nology, Vol.  6,  pp.  279-286,  (1973^. 

7.  Hammond,  D.  C. , Jr.  and  A.  M.  Mellor;  Analytical  CEilculations  for  the 
Performance  and  Pollutant  Emissions  of  Gas  Turbine  Combustors.  Combustion 
Science  and  Technology,  Vol,  4,  pp.  101-112,(1971). 

8.  Roberts,  R.,  L.  D.  Aceto,  R,  Keilback,  D.  P.  Teixeira,  and  J.  M.  Bonneli: 
An  AnalytlcaJ.  Model  for  Nitric  Oxide  Formation  in  a Gas  Turbine  Ccmbtistion 
Chamber.  AIAA  I^per  No.  71-715,  (1971). 
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10).  Chemistry  Is  frequently  modeled  by  acsuminK  equilibriu.m  hydrocarbrr. 
fuel  decomposition  and  two  phase  flow  effects  are  seldom  con.sidered.  In  some- 
more  recent  modeling  attempts,  for  example,  Fletcher  and  He:/wood  (ftef.  5) 
Hajnnond  and  Mellor  (Refs.  6 and  7),  the  stirred  reactor  concept  is  employed 
to  assess  the  effect  of  residence- time  on  combustion  behavior  and  to  predict 
pollutant  emissions  in  gas  turbines.  Droplet  vaporization  and  t'urning  were 
neglected  in  these  studies;  however,  a quasi-global  finite-rate  hydrocarbon 
combustion  mechanism  was  employed  by  Hammond  and  Mellor  to  model  the  chemis- 
try. In  related  work,  Roberts,  et  al. , (Ref.  8),  in  an  attempt  to  predict 
nitrogen  oxide  formation  in  gas  tvirbine  combustors,  subdivided  the  combustor 
into  three  regions:  one  corresponding  to  the  central  recirculation  portion 

of  the  upstream  zone;  a second  representing  the  flow  region  surrounding  the 
recirculation  zone  which  was  interpreted  to  be  a one- dimensional  react- ng  zone 
and  the  third  downstream  zone  modeled  as  a one- dimensional  region.  Both 
finite-rate  sind  equilibrium  hydrocarbon  chemistry  models  were  conside:.  3d.  It 
is  interesting  that  little  difference  in  the  predicted  NO  levels  was  noted  in 
their  results  between  the  equilibrium  and  finite-rate  hydrocarbon  cases.  A 
more  recent  analysis  directed  toward  low  power  application  by  Hosier,  et  al. , 
(Ref.  9)  basically  extended  the  work  of  Ref,  8 through  the  use  of  a more 
sophisticated  finite-rate  hydrocarbon  chemistry  model  obtaining  trends  in 
eigreement  with  experimental  data.  The  modular  approach  proposed  by  Edelmar. 
and  Economos  (Ref.  lO)  is  an  attempt  at  formulating  a general  analytical  pro- 
cedure for  predicting  combiistor  behavior  by  treating  the  various  critical 
combustor  processes  on  an  individual  oasis  or  coupled  as  a function  of  oper- 
ating conditions.  Difficulty  with  the  approach  lies  with  its  method  of 
accounting  for  recirculation  (a  stirred  reactor  is  presently  used)  and  its 
Inability  to  provide  a \inifled  description  of  the  burner  under  a given  set 
of  operating  conditions. 

The  foregoing  methods  are  lacking  primarily  in  their  ability  to  properly 
account  for  mixing  phenomena  occurring  in  the  reverse  flow  recirculation 
regions  of  combiistion  devices.  It  has  recently  become  feasible,  however,  to 
treat  more  rigorously  flows  having  recirculation  zones,  by  n\imerical  solution 


9.  Mosler,  S.  A.,  R.  Roberts,  and  R,  £.  Henderson:  Development  and 

Verification  of  on  Analytical  Model  for  Predicting  Ilnissions  from  G^-s 
Turbine  Engine  Combustors  During  Low  Power  Operation.  Ulst  Meeting  Pro- 
pulsion and  Energetics  Panel  of  AgARD,  (1973). 

10.  Edelman,  R.  and  C.  Economos:  A f-lathematical  Model  for  Jet  Engine 

Combustor  Pollutant  Emissions.  AIAA  I^per  No,  71-711+,  (1971). 
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of  elliptic  equations  governinj.''  combnst.in..;  i’iovs.  r'r.i-  '■xaniple,  one  such 
numerical  method  based  on  an  ext>licis  point -by- poin*  r-- J.Jjjcation  procedure 
has  been  siiggested  by  Gcsaati,  et  al,  , (Ref.  11 ).  .'...bseqve.-itiy , the  Imperial 
College  grouii  has  developed  a line  relaxation  procedure  for  solution  of 
elliptic  partiaO.  differentiea  equations  (Refs,  12  arid  13).  Several  applica- 
tions of  this  procedure  to  combuctiru'  flows  have  been  made  in  recent  years 
with  emphasis  on  the  irdluence  of  both  combustion  chemistry  models  and  tur- 
bulent fluctirations  (Refs.  lU  fuid  15 )• 

For  several  years  Improved  numerical  procedures  have  been  under 
development  at  UTRC  for  solving  com’ou.sting  flows  containing  recirculation 
zones.  One  of  these  IJTRC  proced'rres  is  an  implici'  computational  scheme,  and 
is  novel  In  tiiat  residuals  fii'e  relaxed  3ia;.ultatieously  tlircugiiout  t.he  entire 
flow  field,  rather  tlxan  one  at  a tl-me,  as  in  the  explicit  point  methods,  or  a 
line  at  a time  as  in  the  line  methods.  Under  a joint  /'u7iPL/FAA  contract,  the 
UTRC  Field  Relaxation  i^lliptlc  Procedijre  (FRE?)  based  on  the  rigorous  solution 
of  the  governing  equations  vrtts  tTii-ther  develojx-'d  and  used  ‘-.o  predict  the 
oerformance  aiid  emission  chaxacterlsf  ic.s  of  con-anniilar  and  annular  gas  tur- 
bine combustors  (Ref.  l6).  Further  development  of  the  proced'ure  and  the  phys- 
ical models  is  presently  being  carried  out  under  FPA  Contract  No.  68-02-1873. 
The  UTRC  method  solves  the  axisymmetric  time -averaged  Navier-Ctokes  equations 
Including  the  effects  of  turbulence,  chemistry,  radiation,  and  droplet  vapor- 
ization. 'Ihe  FRl'.P  code  h/ic  given  reasonable  predictions  for  combiostor  flow 
fields  which  nave  axial  sjanmeti'y;  h'cwcver,  significant  circumferential  asymme- 
try Is  present  in  many  aircraft  cornoustor  flow  fields  oiid  a realistic  approach 


11.  Gosma.-i,  A.  P. , W,  M,  Pun,  A.  K.  Riinchal.  D.  P.  fpalding,  at^d  M.  Wolfstein 
Heat  and  Mass  Trarisfer  In  Recirc\J_atir.c  Flews.  Academic  Itress,  New  York, 

(1969). 

12.  Gosman,  A.  D.  and  W.  M.  'I\>n:  Lecture  fiOtes  for  Coui’se  Entitled  "Calcifia- 

tlon  of  Recirculating  Flews",  Report  No.  iiTU/74l2,  Department  of  Mechan- 
ical Engineering,  lmi>erial  College,  London,  (197^). 

13.  Patankar,  S.  V.  eind  D,  B.  Spalding:  A Calculation  Procedure  for  Heat, 

Mass  and  Momentum  Transfer  In  'rtiree-Mmensional  1-a.rabolic  Flows.  Int. 

J.  Heat  Mass  Transfer,  Vol,  15,  p.  1787,  (1972). 

14.  Khalil,  £.  E.,  D,  B.  Spaldlrig  and  J.  H,  NTdielaw:  The  CatLctdation  of 

Local  Flew  Properties  in  IVo-Pijuensional  Furtiaces.  Int.  J.  Heat  Mass 
Transfer,  Vol,  I8,  pp.  775-791,  (1975). 

15.  Gosman,  A.  b. , F.  C.  Lockwood  and  S.  A.  Syed;  Prediction  of  a Horizontal 
Free,  Turbulent  Diffusion  Flame.  li'ocee'lings  Sixteenth  Symposiumi  (inter- 
national) on.  Combustion,  Ihe  Combustion  Institute,  (1976)  to  be  published 

16.  Anasoulis,  R.  F. , H.  McDonald  and  R.  C.  Puggeln:  Development  of  a 

Combustor  Flow  Analysis,  lart  I:  'Ttieorcddoal  Studies.  Air  Force  Aero 

ft’opulslon  Laboratory  Report  AFA1'L-'’'R- 73-98,  f'^ar'  ],  January  I974. 
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au£t  consi-.ler  tliis  corrplication,  ouch  a gcneraJ.  approach  in-/ol.vcu  solution 
of  the  time  averaged  three-dimensional  Navier-Stokf^s  equat'onn  inclurlinf:  tur- 
b\ilence,  chemical  reactions,  radiation  transport  and  droplet  vaporization. 
Recently,  for  instance,  Pataiikar  and  Spalding  (Refs.  17  and  l8)  have  developed 
a line  relaxation  procedure  for  conqmtation  of  steady  ttiree- dimensional  com- 
busting flows  in  cartesian  or  cylindrical  polar  geometries;  however,  this 
procedure  is  not  generally  available  and  few  details  of  this  tectmique  have 
been  disclosed  at  the  present  time. 

Consequently,  under  Piase  I of  the  present  contract  a tiiree-dimensional 
combustor  flow  analysis  progi-am  (Ref.  19)  was  developed  by  extending  an 
existing  and  relatively  efficient  uTRC  tiiree-dimenslonal  Kavier-Stokes  cal- 
culation procedure  (Refs.  20  a:':d  21)  so  that  it  would  be  able  to  compute  com- 
busting flows.  In  particular,  the  procedure  developed  included  a simple 
mixing  length  turbulence  model,  a x’Seudo-kinetic  hydrocarbon  chemistry  model, 
a liquid  droplet  vaporization  model,  a single  freqttency  radiation  model  and 
a fl.nite  rate  nitrogen  oxide  chernistr:/  model. 

The  preliminar;,'  resiolts  obtained  under  Phase  I for  a rectangular 
combustion  chamber  demonstrated  the  feasibility  of  calculating  ttoee- 
dimenslonal  time-dependent  combusting  flows  using  the  MINT  (f-fulti dimensional 
Implicit  Nonlinear  Time-depende.nt ) procedui'e.  Experience  from  the  Phase  I 
study  and  from  the  two-dimensional  EPA  sponsored  FREP  program  pointed  out  the 
severe  shortcomings  of  the  mixing  length  turbulence  model  due  to  the  consider- 
able difficulty  in  specifyi.ng  a reasonable  mixing  length  distribution  even 
for  relatively  simple  combustors.  In  addition  the  problem  of  significant 
truncation  error  because  of  the  relatively  coarse  computational  mesh  (17  x 


17.  I^itankar,  S,  V.  and  D,  B,  Spalding:  A Computer  Model  for  liiree- 

Diraenslonal  Flow  in  IXirnacec.  lYoceedlngs  Fourteenth  Symposium  (inter- 
national) on  Combustion,  Ilie  Corabtustlon  Institute,  pp.  606-6lk,  (1973). 

18.  Patankar,  S.  V,  and  D.  B.  Spalding:  Simultaneous  lire  notions  of  Flow 

fttttern  and  Radiation  for  'Utree-Dimensional  Flames.  Heat  '["ransfer  in 
Flames,  Edited  by  N,  Afgan  and  J.  f^er.  Scrlpta,  Washd.ngton,  (197^). 

19.  Glbeling,  H.  J. , H.  McDonald  and  W.  R.  Eiriley:  Development,  of  a 'I'hree- 

Dimensional  Combustor  Flow  Analysis,  Vol.  I:  Theoretical  Studies. 

Air  Force  Aero- Propuls  ion  Laboratory  Report  AFAPI.-'l’R-75-59,  Vol.  1, 

July  1975. 

20.  Briley,  W.  R.  and  H.  McDonald:  An  Implicit  Numerical  Method  for  the 

Kultidlmenslonal  Compressible  Navier-Stokes  F;quations.  United  Aircraft 
Research  Laboratories  Report  M9II363-6,  November  1973. 

21.  Briley,  W.  R. , H.  Mcltonald  and  H.  J.  Glbeling:  Solution  of  the  M\ilti- 

dlmenslonal  Compressible  Navier-Stokes  Equations  by  a Generalized  Implicit 
Method,  United  Technologies  Research  Center  Report  H75-911363-I5, 

January  1976. 


4 


9 X 17  grid')  which  wan  employed  t'ecam»  apf'.arent  riarin^;  r'f.aue  I of  this  i^rograra. 
In  three-dimensional  con;t)US’'ing  flow  . alc.daM'.ns  at  piv-jen^  the  computer  r’un 
times  are  still  too  long  to  permi’.  }ner<M3ing  Uie  numoer  oi'  mesh  points  sig- 
nificantly to  attain  the  desired  accvu'acy,  Tne  Hiase  I resoits  also  indicated 
the  need  for  a practical  finite-rate  chemistry  model  for  premixed  aircraft- 
fuel  and  air  mixtures  which  could  t>e  incorporated  into  f-Ji  elliptic  computa- 
tional procedure  such  as  MINT,  llie  Latter  problem  Is  one  of  considoraule 
complexity,  and  it  seemed  appropriate  defer  work  in  this  area  until  a 
future  study. 

Therefore,  the  Phase  II  effort  milftr  the  present  contract  was  directed 
toward  incorporation  of  both  a two-equiition  turbul.ence  model  and  certain 
approximate  nonlinear  difference  formulas  into  the  MINT  combustion  code.  Ilie 
so-called  nonlinear  difference  formulas  refiresent  an  ad  hoc  adaptation  of 
five-point  central  difference  formulas  into  the  three-point  implicit  differ- 
ence framework  of  the  MINT  code:.  In  order  to  alleviate  the  problem  of  trunca- 
tion eirror  associated  with  coarse  corqiULational  grids.  Wlale  the  varicAis 
formulas  considered  gave  acceptable  results  in  noru'eactlng  flov/s,  unfortu- 
nately no  convergent  solutions  v;ir,h  the  nonlinear  i’orr.'iulas  were  obtained  for 
combusting  flows . Therefore,  the  namericai,  predictions  presented  herein  were 
obtained  using  conventional  central  difl'erence  formulas.  As  a consequence  of 
being  ’uriable  to  refine  the  computations,  grid  due  tr  coir.puter  run  time  limi- 
tations, the  resulting  combusting  flow  caiculations  witl.  conventional  central 
difference  formulas  almost  certainly  contain  significant  truncation  errors, 
thus  comparisons  with  experiment  can  only  be  rtwed  qualitatively.  In 
addition,  tne  inlet  conditions  for  the  eXierir.eiK.  are  knovnr.  oi.Iy  approximately, 
further  confusing  the  coiuparisori  lenwe'-n  predictior.r;.  and  measurement.  On  the 
positive  side,  it  was  possible  tc>  perfv>rm  stal’le  ca;  ■■ilia’. ions  for  a very  com- 
plex reacting  flow  using  a two-equatior.  isodei  of  tiu-bulence . Ttie  principle 
obstacle  to  a thorou^i  evaliaation  of  the  procedure  is  simply  insufficient 
grid  points  together  with  insufficient  data  on  inlet  con  ILtions.  Tlie 
increased  storage  requirements  of  tlie  additional  grid  points  required  to 
define  the  flow  field  accurately  do  not  pose  a parti ciali.u’  problem  in  view 
of  the  efficient  use  of  mass  storage  possible  with  the  MINT  code. 
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'•  MOI-fL 


T.rr. T ; Eq'.iatiO.iE 

The  flow  in  combiis'i  . : J \ ' V:.ov^'  '■-o  be  predon;'. rai.!  ,■  I 'U'b"Ie.Tit. . To 

account  for  this  turbiile'. : r '.-i-n'rr  ‘/-.c  solution  oi'  the  tic.t -hvc  ra>/.e 

iN'avier-Stckes  equations,  a : r,  h'::..'-.-  is  introduced  hefine  an  effective 

viscosity.  A review  of  tua'sulence  sicdcis  is  available  in  the  literature  (e.g. 
Refs.  22  to  2U),  Prar.dtl  (Ref.  25.  was  perhaps  the  first  to  inti'oduce  a turbu- 
lence model  when  he  postulated  that  *:ht-  t Ime-averap.ed  shear  ■.  tre?s  and  the 
time-averaged  velocity  gradient,  ai'e  proportional  as  in  lar/iir.ar  flow,  and  that 
tne  length  scale  (the  so-called  mixing  length)  vshich  enters  the  relationship 
is  proportional  to  the  ‘urbuir.s  shcai’  re,,uon  thickness.  Prandf  ‘ mixing 
length  model  has  oeen  empic;y>  s --;ac"Or.of.:.l  Ly  by  a number  of  ir  .-e-sti gators 
(e.g.,  Refs.  26  to  3'u;  in  a ri ' t..,  f problems-  primarily  ir-volving  turbulent 
fi.aw  along  '«/alls  and  in  free  rui-hu*e.it  flo-ws.  A disadvantage  of  t;ie  mixing 
length  model  is  that  it  is  an  eqall.Isriu;;.  model  (;.e.,  turbulence  is 
assumed  to  be  produced  and  dissipated  locally)  and  it  requiies  an  ad  hoc 


22.  l.aunder,  E.  E.  and  D.  B.  Spalding;  Mathematical  Models  of  Turbcilence. 
Academic  Press,  London,  1972, 

23.  Harlow,  H,,  ed.;  Turbul enc-»  Trar.snort  Modeling.  ATAA  Selected 

Reprint  Series,  Vol.  I07". 

24.  Launder,  B.  E. ; Itrogre:'  s fr  the  Modeling  of  T-urbulenf  '"i-aricport. 

Notes  for  Coturse  on  Turbulent  Reci.rculatJ  ng  Flows -T'redi  'ti on  and 
Measurement,  Pennsylv'ani State  University,  July  197o. 

25.  Prandtl,  L.:  Bericht  Uber  Urjtersuchingen  7,-ir  Aus,.:ebildeten  IPirbulen/:. 

zm-:,  Vol.  5,  1925,  p.  136. 

26.  Patankar,  3.  V.  and  P.  Epnldir./y  Heat  and  Mu;,s  Trat.r.fer  in  Boundary 
Layer.s.  Intertext  Books,  Ixiraon,  1970. 

27.  Malse,  G.  and  H.  McDonald;  Mi  yin.'-  i.engbh  and  Kinematic  Eddv  Viscosity 
in  a Compressible  Boundar  Lay'.r.  /MAol  Journal,  Vol.  ’yoP,  pp.  73-SO. 

28.  McDonald,  H.  and  T'.  J,  Camarata;  An  Extended  Mixing  Length  Approach 
for  Computing  the  Turbulent  Bc-’ondary  I.ayer  Development.  Proceedings 
of  the  AF03R-IFP-Stanfo"cd  Conference  on  Boundary  Layer  Prediction, 

1968. 

29.  Williaimson,  J.  W.:  Ati  Extension  of  Prandtl's  Mixing  Length  Tiieory. 

Applied  Mechanics  and  Fliiids  Engineering  Conference,  A3tdE,  June  l(l69. 

30.  Lllley,  D.  0.:  Prediction  of  Inert  Turbulent  Swirl  Flows.  AlAA  Paper 

No.  72-699,  1972. 


6 


mixing  length  distribution.  Some  ol‘  the  r.hortcoiningo  o;'  the  mixing  length 
model  have  been  overcome  for  many  cases  of  Interest  by  tlie  introduction  of 
various  multiequation  transport  models  ol'  turbulence. 

Many  of  the  two-equation  turbulence  models  employ  the  Prandtl-Kolmogorov 
formula  for  specification  of  the  tui’uulent  viscosity, 


tl 

Re' 


(1) 


where  k is  the  turbulence  kinetic  enerfr/  and  Z iB  a length  scale  of  the 
t\irbulence.  'Phis  relation  follov.'s  from  dimensional  argisTients  for  t'urbul  ent 
flow  described  by  the  two  parameters,  k and  1.  A major  advantage  of  a two- 
equation  turbulence  model  compared  to  a mixing  lengnh  model  is  the  fact  that 
the  length  scale  (as  well  as  the  turbulence  kinetic  energy)  is  determined 
from  transport  equations,  whereas  in  irdxlng  length  model.s  the  length  scale 
is  determined  from  an  ad  hoc  algebraic  expression.  Successful  use  of  the 
mixing  length  model  requires  an  a pi-ioxu  specification  of  the  turbulence 
length  scale.  While  a realistic  assun.ption  can  be  made  for  certain  flov,’s 
such  as  shear  layers,  in  many  flows  of  interest  (such  as  interrial  recircu- 
lating flows)  the  choice  of  a proper  turbulence  length  scale  is  not  obvlo’us . 
Since  the  turbulence  length  scale  emerges  from  the  solution  in  a two-equation 
model,  these  models  are  more  likely  to  give  accurate  prediction.",  over  a 
wide  range  of  geometric  and  flow  conditions  with  the  same  empirical  constants. 

It  shoui-d  be  noted  that  tiie  two  equation  models  employ  tire  eddy-viscos ity 
formulation  for  the  Reynolds  stresses  as  in  the  mixing  length  model,  i.e., 


P 


Re 


(2) 


Hence,  this  formulation  still  suffers  from  the  physical  siiortcoming  that 
there  is  zero  Reynolds  stress  wherever  the  velocity  gradient  is  zero.  In 
addition,  the  eddy  viscosity  form'ilatlon  is  isotropic  which  may  be  incorrect 
in  many  three-dimensional  and  swirling  flows.  However,  for  practical  calcu- 
lations of  turbulent  internal  recirculating  flows  there  are  no  other  available 
transport  niodels  which  areas  suitable  or  even  as  relatively  well  developed. 

In  the  present  MINT  computer  code  there  are  three  mixing  Length  models  and 
one  two-equation  transport  model  of  tiurbulence.  Each  mixing  length  model 
can  be  used  individually  or  criteria  can  be  specified  such  that  the  appropri- 
ate mixing  length  model  will  be  utilized  in  various  regions  of  the  flow  field. 
These  models  will  now  be  discussed  in  some  detail. 
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The  mixing  length  turbuJenct'  mu-i' J:  -‘.inployod  in  Uw;.  ana  I vf:  is  are  bared 
on  mixing  length  dintributionr  oUgf-n  t\v  Willi, nmron  (b,  '.  , Id  Iley 

(Ref.  30)  , and  van  Driert.  (ib-t  . J ; . Wi  Iliitiitiion'f,'.  modoi  ? n'iicable;  to 
to  flows  in  channels  and  piper.  i/.Ilcy'.;  model  is  appiicab:-.-  to  uriconfined 
jet  flows,  and  the  van  Driest  model  is  valid  for  flows,  near  solid  wells. 

The  mathematical  fonii  of  the  expression  for  the  tuj’buleni  viscosity  follows 
from  Ref.  32: 


Re 

where  e is  the  mean  flow  rate  of  strain  tensor 

e-  1/^  [(Vu)  ->  (Vu  f] 

and  the  mixing  length  £ is  givt'n  >_y  either  Wi  lliarrr;on ' s 


(3) 


or  Lilley's  model 


7-=oi4(p-lexo(i--) 


^ = 006(^,1 


O w 


= 0 05 


or  van  Driest ’s  model 


£ : #C  y 


-exp(-y7A*)] 


(5) 

(6) 


(7) 


In  the  above  expressions  for  £,  r^  is  the  distance  from  the  wall  to  the 
centerline  of  the  duct  or  pipe,  y is  the  distance  from  the  point  in  question 
to  the  nearest  wall,  4 is  the  half-width  of  a shear-layer  like  region  as 
defined  by  Eq.  (6)  and  k and  are  the  von  Karman  constant  and  van  Driest 
damping  coefficient,  respectively.  The  constants  in  Eqc.  (‘J-7)  have  been 
determined  by  comparison  of  theory  and  experimental  data  for  the  class  of 
flows  under  consideration. 


31.  van  Driest,  E.  R.:  On  Turbulent  Eiov,'  Near  a Wall.  Jouriiul  of  the 

Aeronautical  Sciences,  November  i '7'  . 

32.  Beer,  J.  M.  and  N.  A.  Chigler:  Cuinburtlon  Aerodynamics.  Wiley,  New 

York,  1972. 


VarioUE  formr.  of  the  two  equation  mo;k:i  of  turbiJ.enco  have  been 
proposed  since  Koinogorov  (Ref.  33)  first  introduced  the  concept  in  1942. 
Most  investigators  have  chosen  the  kinetic  of  tun/alence,  k,  as  their 

first  variable.  However,  there  is  a wide  diversity  of  choice  as  to  the 
second  variable  to  be  used.  In  general,  eacli  investigator  chose  a second 
variable  which  he  felt  was  appropriate  to  the  physical  description  of  tur- 
bulence. For  instance,  Kolmogorov  chose  as  Ids  second  variable  a quantity 
which  was  proportional  to  the  mean  frequency  of  the  most  energetic  motions, 
while  Spalding  (Ref.  34)  and  Safiir.an  (Ref.  35)  chose  a quantity  that  repre- 
sented the  time-average  square  of  the  vortlcity  fl.uctu.ations . Another 
commonly  chosen  second  variable  has  been  the  turbulence  kinetic  ener(:;y 
dissipation  rate,  e,  which  was  selected  for  this  investigation.  An  advan- 
tage of  using  the  dissipation  rate  of  tiurbiolence  kinetic  energy  for  tiie 
second  equation  is  that  the  dissipation  rate  appears  directly  in  the  turbu- 
lence kinetic  energy  equation,  and  tiiat  an  equation  for  ^ can  be  readily 
developed.  Since  derivation  of  the  equations  for  k and  ^ are  lengthy  and 
have  previoiosly  been  presented  in  the  open  literature,  there  derivations 
are  not  repeated  in  the  present  report.  Ihe  appropriate  transport  equations 
for  turbulence  kinetic  energy  and  energy  dissipation  rate  valid  at  high 
Reynolds  numbers  are  taken  directly  from  Launder  and  Spalding  (Ref.  36). 


The  turbulence  kinetic  eneiq^yr  equation  in  vector  notation  is 


d(Pk) 

dt 


(8) 


where  the  first  two  terms  on  the  right  hand  side  represent  convec-ion  and 
diffusion  of  turbulence  kinetic  energy,  respectively;  the  thira  and  foiuth 
terms  represent  the  generation  (due  to  shear  forces)  and  dissipation  of 
turbulence  kinetic  energy,  respectively.  The  equation  for  the  dissipation 
rate  of  turbulence  kinetic  energy  is 


(31 


■ - V (/>S.l  + ^ 7.)+ C,  ^ f- ( 2l  l ) 


(9) 


« 


33.  Kolmogorov,  A.  N.;  Equations  of  Turbulent  Motion  of  an  Incompressible 
Turbulent  Fluid.  IZC.  Adak.  Naut.  SSR  Ser.  lliys . VI,  No.  1-2,  56, 

1942. 

34.  Spalding,  D.  B.:  Tie  Prediction  of  Two-dimensional,  Steady  Turbulent 

Flows.  Imperial  College,  Heat  Transfer  Section  Report  EIVtn/a/16,  1969. 

35.  Saffman,  P.  G.:  A Model  for  InJiomogonous  Turbulent  Flow.  Proc.  Roy. 

Soc.  Ser.  A317,  1970,  p.  417. 

36.  Launder,  B.  E.  and  D.  B.  Spalding:  Tiie  Numerical  Computation  of 

Turbulent  Flows.  Computer  Methods  in  Applied  Mechanics  and  Engineering, 

Vol.  3,  1974,  p.  269. 
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where  the  function  of  the  terms  is  analaf.'ou3  to  that  in  the  t'orbulence 
kinetic  energy  equation.  In  Eqs.  (8)  and  (9).  e is  the  mean  flow  rate  of 
strain  tensor  defined  by  Eq.  (4).  Using  dimensional  arfaiments  the  I'randt] - 
Kolmogorov  formula,  Eq.  (l)  may  be  written  as 

Re  V f 


which  implies  a turbulence  length  scale  or  ''mixing  length"  nay  be  defined 
as 


(11) 


The  constants  appearing  in  Eqs.  (8-10)  have  been  evaluated  by  Jones 
and  Launder  (Ref.  37)  suid  those  values  will  be  used  in  this  work. 

= 0 09 

C,  : I 55 

C2=20  (12) 
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Boundary  Conditions 

Fomulation  of  suitable  bo\indary  conditions  to  be  employed  with  the 
k-<?  turbulence  model  has  proven  to  be  critical  in  recirculating  flows. 
Launder  and  Spalding  (Ref.  36)  advocated  setting  the  value  of  dissipation 
rate  (e)  at  one  grid  point  away  from  the  wall  consistent  with  an  analytic 
wall  function  formulation.  Tiie  appropriate  relation  is  obtained  by  con- 
sidering the  turbulence  length  scale  in  the  near  wall  region. 

The  eddy  viscosity  mixing  length  model  for  is 


37.  Jones,  W.  P.  and  B.  E.  Launder:  The  Prediction  of  Laminarization  with 

a TVro-equation  Model  of  Turbulence.  Int.  J.  Heat  Mass  Transfer,  Vol.  15, 

1972,  p.  301. 
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rare* 
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(13) 


n : 1/2 

— =/)i2(2e  e) 

which  is  vBlid  near  a wali,.  Tl'ie  constaxit  C^,  Eq.  (11),  may  te  evaluated  by 
assuming  turbulence  equilibrium,  i.e., 


^ - - 2 - -i/Z 

€-  (2e  e)  = i (2e  e) 


Using  Eq.  (U)  and  the  Prandtl-Kolmogorov  fonnula  yields 


i,'/2 


Ri'-/’^^2e§)  =c^- 


(14) 


(15) 


so  that  one  finally  obtains 


Cf-  C|i 


3/4 


(16) 


Then  Eq.  (U)  provides  a relation  for  r near  a wall  where  the  mixing  length 
£ Is  assumed  to  be  known: 


« =C, 


3/4  kV2 


(17) 


The  value  of  kinetic  energy,  k,  at-  the  first  grid  point  away  I'rom  the  wall 
is  obtained  by  solution  of  the  k-transport  equation,  along  with  the  approxi- 
mate condition  ^^k/^^n  = 0 applied  at  that  'point. 

Results  from  a two-dimensional  study  being  performed  under  ERA  sponsor- 
ship at  UTRC  have  led  to  a revision  in  the  wall  function  concepts  described 
in  Ref.  I9.  The  modified  wall  function  approach  appears  to  be  less  restric- 
tive than  the  method  described  in  Ref.  19,  ^nd  seems  to  give  good  results 
for  a variety  of  flow  configurations.  Tl.e  basic  premise  is  that  at  one  grid 
point  away  from  the  wall  the  normal  derivatives  of  turbulent  viscosity  and 
velocity  should  be  consistent  with  a wall  fur.ction  formulation.  It  can  be 
shown  from  the  analysis  of  van  Driest  (Ref.  "^-l)  that  the  general  relations 
which  are  valid  in  the  near  wall  region  are 


4-  •♦'i  4*. 

u - u ( y ) 


(18) 


_L  iiU  . I 

mt  ay  ' Q(y.y^) 


, au  ' du*’  ^ 


(19) 

(20) 
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where  u i2  the  total  velocity  conipoaf'r.t  parallel  to  the  wal  ! , y,,.,  ir;  the  t .ri 
^ent  viccooity,  and  y is  the  dintarice  i'ror.  the  wall.  Also,  y^  \ar.d  ar' 
defined  in  the  conventional  imrinor  a" 


(:  J j 


and  tiie  friction  velocity  u i:-  yiveri  t,y 


The  f’jnction  c(y,  y'*')  in  Eq.  (I'j)  varier  from  y/P  in  the  lan.lnar  su:layei- 
to  y in  the  logarifnmic  region.  Ti.e  i’unctional  fonri  for  u"^  Implied  ;y 
Eq.  (l8)  coaid  be  obtained,  for  example,  from,  the  anaiysic  of  van  Driest 
(Ref.  31}  or  from  the  more  empirical  relation.c  given  by  VAalz  (Ref.  38). 

In  the  .-0.^ari.nmLic  velocity  profile  region  (y"*  ^ 80),  the  conventional  law 
of  the  wall  can  be  written 


i 


ar:d  relations  (ih-RO;  redace  to 

— r' . ' _i_ 

ay  ■ dy‘  ^ ay2  ' y 

These  relations  (25)  seir/e  to  specify  a wall  value  of  turbulent 
'^^“-t'-ity  and  a wal  i slip)'  velocity  such  that  the  normal  derivatives  are 
consistent  with  the  logarithmdc  wall  function  formulation.  Note  that  Eq. 
(2>/  IS  independent  of  the  constants  cp  and  C2  appearing  in  Eq.  (2k),  arid 
hence  this  approach  appears  to  be  less  restrictive  than  using  Eq.  (PU) 
directly.  Tl.ese  relations  have  been  employed  in  the  calculations  carried 
out  'under  idiase  II  of  the  present  contract  without  difficulty,  although  flow 
resol jtion  in  the  region  very  near  the  wall  is  sacrificed  to  attain  accuracy 
in  the  central  region  ol'  the  flow  field  where  the  comtiu.sf.i  on  firfir-e.sse.s  are 


38.  Walz,  A.;  Boundary  Layer:;  of  Flow  and  Terrifie rat.uri ■ . M.I.T.  I’re::;-, 
Camt.rid,''e,  Massachusetts , p.  115>  I'X’A'. 
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co:. . hov.’ever,  1:’  accirato  caLculationc  in  th^j  vail  region  arc- 
required  ai  a *atcr  date,  reiinciner.tr  to  tl;e  valL  furiCtion  anproach  ina.v 
L€  in[  i-ernented  earily  in  t.nt  precent  eonputati ona L proced  jre. 


Ariotiier  problem  arlolri''  in  t.he  ni-actical  application  ol’  a two  equation 
^ urnulcnce  iriodel  which  han  not  beeri  widely  di  nc'tnaed  i.-:  the  difficulty  of 
opecifyin,;  suitable  function  boundary  coriditionn  for  k and  «:  at  a computa- 
tiot.al  inlet  yjlane.  In  modelinr  a general  exrjcrimentai  conf iruration , the 
pr'yi.ler.  ol‘  cpecifying  k and  r at  the  inlet  plane  in  complicated  by  the  lack 
of  detailed  inlet  velocity  and  turbulence  profiles . iience,  it  is  necessary 
tc  a-ssur-.e  the  inlet  velocity  profile  and  then  construct  consistent  profile.n 
for  k and  e therefrom.  I.n  general  one  possibility  is  to  assume  turbulence 
equilibri’um.  and  a mixing  length  distribution,  i.e.,  Eq.  (l4).  Then  usit.m 
Eq . IC,  , 


f^(2e  e) 


(2e) 


arid  e foliovrs  frorr.  Eq.  fl^t).  A djffl--.!ty  v/ltb  this  ar. preach  occiurn  wlieii 
the-  Vf-lo-ity  firofiic-  beco.m.er.  un:  forr-i,  the  j'atf-  o!’  .".trairi  teri::or  fi  ' 

ar  proache;,  zero.  It  i.-:  kriovrn  ttiut  a r, t ri-rirri"  tur.UuJe-nce  level  exl.nt: 

i virtually  all  uru  1‘orr;  flow;;,  ::  i;.  t i r.  treat  a ::ri<.-ci  fi  e-u  fr-''"-  .-trerim 
turr.iilence  kinetic  e-n'-rgy,  Xp..  ne  au-;»--J  to  E'i.  ffh)  to  r»roviue  tlio  total 
in*et  turn  iience  level, 


k = kp;  ♦ 


<^(2l  I) 


(2Y) 


where  / is  still  the  snecifieu  .m.lxit..-  len-'th.  Tr.c  addition  of  free  stream; 
turbulence  implies  that  equilibrium,  no  longer  exists  in  ger.eral,  and  e 
follows  from  Eq.  (17).  The  addition  of  free  stream;  turh  iience  as  indicated 
in  the  above  procedure  .scemr.  to  rirovide  reasonable  values  for  the  turbulence 
variables  k,  <'  and  and  the  resiiiting  .-olutionr  exdiibit  qualitatively 
resonablc  behavior. 


Numerical  gori::  iU(  rati  or.,-; 


Efficin-rit  ;:fil  litl  on  nl'  tbr-  k--  mo'i<  I gov-rnirig  <-qunt  i or... , i/p  . 

Vig/rthf-r  '*f  1 tb  tb<-  ''oii;;<  ■ rvati  fin  f?quat  i 'jn:;  '.f  , rTKmi'-rd.ujM,  '■tji-r'i'v  ati-l 

Specie.';  desf-r  1 r.f.-'l  in  P'.-f.  11  rc-qui r'.-;;  /)i-oi<er  treutmerit  'jf  I, he-  .'-.ourci.,-  tf.-nm; 
in  Eq::.  ('jj  arid  fl).  'llie  calculation  priicedorf;  of  Hef.  1''  i-mjil(iycd  a mixin. 


13 


length  turbulence  model  and  an  eddy  viscosity  evaluated  from  the  knovm  solu- 
tion at  time  level  t^.  This  procedure  has  also  been  'osed  with  the  turbu- 
lence model,  Eqs . (8-10).  Tlie  source  terras  in  Eqs.  (8-9)  may  then  be  linear- 
ized in  accordance  with  the  technique  discussed  in  Refs.  (19  to  21).  Of 
course,  it  is  then  necessary  to  solve  the  k-€  equations  in  a coupled  manner. 
However,  since  k and  e serve  only  to  specify  the  turbulent  viscosity  through 
Eq.  (10),  it  is  in  fact  possible  to  solve  the  k-e  finite  difference  equations 
after  determination  of  the  flow  field  for  a given  time  step.  This  approach 
requires  considerably  less  computer  run  time  than  the  alternate  procedure 
of  solving  the  complete  system  of  algebraic  equations  in  a coupled  fashion, 
and  has  been  employed  successfully  in  computing  the  results  contained  in 
this  report. 


Ik 


SECTION  III 


SPATIAL  Om-TIRENCE  FOKWJU'TION 


In  computing  solutions  for  high  Reynolds  numbers  with  a conventional 
centered  second  order  accurate  spatial  derivative,  it  is  often  found 
necessary  to  add  a form  of  artificial  viscosity  or  dissipation  for  the 
axial  flow  direction.  Artificial  dissipation  in  some  form  is  often  useful 
in  practical  calculations  to  stabilize  the  overall  method  when  boundai-y 
conditions  are  treated  inaccurately,  when  coarse  mesh  spacing  is  used,  or 
in  the  presence  of  discontinuities.  The  need  for  artificial  dissipation 
artsee  in  certain  instances  when  centered  spatial  difference  approximations 
are  used  for  first  derivative  terms . The  use  of  arid  ficial  dissipation  is 
thus  a matter  of  spatial  differencing  technique,  and  is  commonly  employed 
either  explicitly  or  Implicitly  in  both  explicit  and  implicit  difference 
schemes.  The  main  problem  with  artificial  dissipation  encountered  in  the 
present  work  is  the  resultant  smearing,  of  the  solution  which  occurs  in  a 
high  Reynolds  number  flow  with  inadequate  mesh  resolution.  The  particular 
form  of  artificial  diffusion  described  here  is  considered  provisional, 
since  the  formal  accuracy  is  lowered  to  first  order  for  the  axial  (x^) 
coordinate  direction.  Another  combination  of  spatial  difference  formulas 
and/or  artificial  viscosity  terms  may  eventmlly  prove  superior.  The 
objective  of  this  phase  of  the  present  work  was  to  investigate  and  utilize 
a spatial  difference  formulation  which  would  yield  stable  solutions  with 
reduced  numerical  smearing.  A preliminary  evalmtion  of  an  alternative 
first  derivative  formula  which  had  some  potential  to  eliminate  the 
artificial  viscosity  was  made,  and  this  alternate  difference  scheme  will 
be  discussed  subsequently. 

The  dissipation  term  used  here  for  the  conventional  second  order 
accurate  first-derivative  is  based  on  an  observation  (Roach,  Ref.  39» 
p.  162)  that  for  a linear  model  problem  representing  a one-dimensional 
balance  of  convection  and  diffusion  terms. 


u,  ^ - -L 
* dx,  Re  dx  2 


(28) 


39-  Roache,  P.  J. : Computational  Fluid  Dynamics.  Herraosa  Publishers, 

Albuquerque,  1972. 
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it  can  be  shown  by  a comparison  with  an  exact  solution  oi'  Eq,  (26), 
solutions  detained  usin</  central  differences  for  the  convection  tern; 
are  well  behaved  provided  the  mesh  Reynolds  number  ju^|Ax^He  is  ^2, 

but  that  qualitative  inaccuracies  (associated  wiuh  boun^ry  conditions)  may 
occur  when  Re^x^  > This  suggests  the  use  of  an  artificial  viscosity  to 
ensure  that  the  local  effective  mesh  Reynolds  number  is  no  greater  than 
two.  Thus,  Eq.  (28)  is  replaced  by 


where 


if  Re . >2 

Ax, 

if  Re.  < 2 
AX, 


(30) 


From  Eq.  (30),  it  is  apparent  that  the  artificial  viscosity  can  be 
made  to  vanish  by  refining  the  mesh.  Since  the  artificial  viscosity  is 
proportional  to  dxq,  solutions  will  have  first-order  formal  accuracy  if 
^ ^ second-order  accioracy  if  Re^xi  ^ 2.  Strictly  speaking, 
the  overall  method  is  second-order  accurate  since  Re^xj  "*  0 as  the  mesh 
is  refined.  It  should  be  remembered,  however,  that  such  asymptotic 
truncation  error  estimates  are  meaningful  only  for  sufficiently  small  mesh 
size;  whereas,  in  practical  calculations  of  complex  flows,  mesh  resolution 
capabilities  have  often,  out  of  necessity,  been  strained.  One  virtue  of 
the  present  formulation  is  that  by  isolating  the  artificial  viscosity  terms 
for  comparison  with  other  terms  in  the  equations,  a nonrigorous  but 
plausible  a posteriori  indication  of  the  first  order  truncation  error  in 
a computed  solution  is  available.  It  is,  of  course,  obvious  that  Eq.  (29), 
upon  which  the  artificial  viscosity  is  based,  represents  a gross  simplifi- 
cation of  the  Navier-Stokes  equations,  and  it  is  primarily  for  this  reason 
that  the  present  formulation  of  artificial  viscosity  terms  is  considered 
provisional. 


Particularly  in  the  high  Reynolds  number  type  of  three  dimensional  flows 
considered  in  the  present  work  the  number  of  grid  points  required  to  keep 
the  cell  Reynolds  number  less  than  two  (Re,;^xi  ^ 2)  becomes  prohibitive  if 
less  than  raultihour  computer  run  times  on  a CDC  6600  are  to  be  avoided. 

In  many  circumstances  of  interest  in  con±iusting  flows,  such  as  diffusion 
flames , the  flow  gradients  and  resulting  physical  diffusion  are  such  that 
the  use  of  an  artificial  viscosity  or  diffusion  coefficient  applied  to  such 
large  gradients  can  seriously  degrade  the  solution  by  the  introduction  of 
significant  artificial  numerical  diffusion.  Various  strategies  have  been 


r 


and  are  being  evolved  to  circumvent  tnie  dil'ficulty.  Of  particular  note  is 
the  revived  interest  (Orszag  & Israeli,  Ref.  Uo)  in  fourth-order  compact 
^ ' difference  formulas  involving  only  ti.ree  grid  points.  Mitchell  (Ref.  Ul ) 

for  Instance,  cites  some  early  applications  of  tlds  technique  within  the 
ATI  framework.  Unfortunately  these  compact  difference  formulas  lose  a 
great  deal  of  their  attractiveness  when  both  first  and  second  derivatives 
in  the  same  coordinate  direction  are  present,  since  it  turns  out  that  the 
first  and  second  derivatives  each  require  a separate  computational  level 
within  the  ADI  procedure.  In  terms  of  computational  effort,  Briley, 

McDonald  and  Gibeling  (Ref.  2l)  estimated  that  the  fourth  order  compact 
difference  schemes  required  about  the  same  amount  of  work  as  did  the  more 
conventional  five  point  centered  fourth  order  accurate  difference  formiulas . 
The  key  point  remains,  however,  that  if  Orszag  and  Israeli's  (Ref.  Uo) 
estimate  is  valid  of  the  possible  reduction  in  the  number  of  grid  points 
for  the  fourth  order  schemes  relative  to  the  second  order  schemes  to 
achieve  comparable  accuracy,  very  significant  benefits  are  to  be  obtained 
from  using  higher  order  accurate  schemes  in  nrultldimensional  problems.  It 
was  not  possible  within  the  constraints  of  the  present  study  to  accon5>lish 
the  necessary  code  modifications  required  to  permit  the  enlargement  of  the 
finite  difference  molecule  in  the  MINT  code  from  three  to  five  points. 

I Nor  was  it  possible  to  introduce  the  conventional  compact  three-point  fourth 

order  scheme  without  a significant  effort  in  view  of  the  additional  levels 
required  within  the  ADI  framework.  Instead  a rather  simple  compact  scheme 
which  could  be  easily  implemented  into  the  present  scheme  was  evaluated. 

This  simple  compact  scheme  will  now  be  described. 

The  centered  fourth  order  accurate  five-point  first  derivative  formula 
can  be  written  for  equally  spaced  mesh  as 

~ 

Next  it  is  assumed  that  between  grid  points  the  function  0 can  be  represented 
I by  some  exponential.  It  then  follows  that 

, /</>.  . = /<^,  (32) 


Uo.  Orszag,  S.  A.  and  M.  Israeli:  Numerical  Sinmlation  of  Viscous  Incompres- 

sible Flows.  Annual  Reviews  in  Fluid  Mechanics,  Vol.  6,  197U,  p.  28l. 

Ul.  Mitchell,  A.  R. : Con^jutatlonal  Methods  in  Partial  Differential  Equations, 

Wiley,  New  York,  I969. 
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thus  allowing  the  five  point  difference  molecule  to  be  compressed  to  three, 
which  after  a little  mar.ip'.ilation  gives 


^ ^ r 2 V '^1-1  ] (33) 

dx  2Ax'  ' l'"  J 


= 8 <^,[i-Ax^S^qb,/6d>,] 


(3U) 


Now  it  is  easily  verified  that  on  linear  combinations  of  sine,  cosine  and 
exponential  functions  the  given  formula  is  fourth  order  accurate.  It  is  also 
apparent  that  the  forrmala  is  nonlinear  and  must  encoiunter  difficulties  as  0^ 
passes  through  zero.  Further  althoi;gh  the  fonniuLa  as  given  is  not  indepen- 
dent of  a translation  by  a constant  in  the  dependent  variable  0^ , translation 
by  an  additive  constant  'Tould  eliminate  the  zero  denominator  problem.  If 
this  eulditive  constant  were  infinite  the  conventional  three  point  second 
order  central  difference  forsiula  for  the  first  derivative  is  recovered,  and 
as  mentioned  earlier  If  t.ne  additive  constant  were  zero  a fourth  order 
accurate  formula  for  sines,  cosineB  a.nd  e.xponentlals,  not  suffering  a cell 
Reynolds  number  of  two  limitation,  is  recovered.  Under  this  effort  an  ad  hoc 
limit  was  placed  upon  the  c or  recti- -.i  ■*.erm  Ax^6^0j, /60j^ ; this  term  was  limited 
to  a value  cf  0.5.  Tlds  procedure  avoids  the  difficult5’-  associated  with  a 
zero  denominator.  Insofar  as  the  nonlinear  aspect  of  the  formulation  is 
concerned  this  was  readily  treated  using  the  time  linearization  concept 
(Refs.  19  to  21)  previously  introduced  to  treat  the  nonlinear  dependent 
variable  terms  In  the  governing  equations.  For  ease  of  implementation  the 
dencminator  term  was  not  linearized  but  treated  explicitly.  This  was  felt 
reasonable  for  a preliminary  evaluation,  since  if  warranted  by  the  results 
this  explicit  term  could  be  treated  more  precisely  at  the  expense  of  some 
tedium. 

Under  the  present  effoii.  the  nonlinear  difference  formula  was  applied 
to  the  continuity  equation  for  a simple  axisynmetric  diffusion  flajae  in  a 
cylindrical  pipe  in  which  there  was  no  reverse  flow.  The  result  of  this 
application  was  a reduction  of  the  mass  flux  error  at  the  compu‘''ational  exit 
plane  from  about  thirteen  (I3)  percent  to  less  ttian  one  percent.  Finally  the 
scheme  was  applied  to  the  three-dimensional  rectangular  combustor  discussed 
In  Section  IV,  but  a convergent  solution  could  not  be  obtained.  The  source 
of  the  numerical  difficulty  was  felt  to  be  related  to  the  still  too  coarse 
mesh  and  to  the  form  of  the  first-derivative  correction  term  In  Eq.  (3^), 
i.e.,  (Ax)^  The  correction  term  can  give  a major  contribution 

to  the  derivative  in  regions  with  large  second  derivatives  of  the  dependent 
variables,  particularly  when  the  variable  0^  is  near  zero,  and  this  may  lead 
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to  erroneous  results.  Finally,  the  one-sided  continuity  equation  wail 
boundary  condition  is  difficult  to  formulate  consistent  with  the  nonlinear 
scheme  Eq.  (3^)»  and  this  also  may  have  contributed  to  tne  numerical  problem. 
Within  the  constraints  of  the  present  program  these  difficulties  could  not 
be  resolved.  At  this  point  the  origin  of  the  difficulties  iias  not  been 
established  precisely  so  it  is  not  even  clear  if  the  above  scheme  could  be 
made  to  work  satlsfactortly  for  reacting  flows  with  recirculation.  Certainly 
the  motivation  to  continue  efforts  to  Increase  the  accuracy  and  minimize  the 
number  of  grtd  points  for  three-dimensional  combusting  flow  calculations 
is  still  present. 
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To  evaluate  the  three-dimenc !■  Mf-. . c inp „•  at  ior.al  pr  'Ce!;  .re  described  in 
the  previous  sections,  ooi.iparic  :i:  wer(;  tr,ad<-  :e*w<:-en  numerical  predictions 
C'btained  from  the  MlNi'  •.•.•de  anu  i!.p  li'':;od  e.’tp':  rLme;.- m uata  taker;  ir.  a rec- 
tangular research  ccmr  is'-  r ty  \]y,'  ■ ■ V,'h’T.r-.ey  ‘'r-.-Taf'  :P?^V.'A)  iornhusti  t; 

Group.  The  P&WA  resul*  .>  c>;;si.s'ea  f • ■■mpera'.  .re  ;nea;'...rexen’ s '^aken  witr.  a 
shielded  thermocouple  pr.,.e  ana  e::,i;u*  ..  measuremerit.'  ; .xr.curned  hyorocart  ns, 
nitric  oxide,  carbon  m.noxide.  carton  dioxide'  take.-,  using  gas  sampling  probes. 
No  other  flow  measurementr  were  avail  a:  i-  f -r  c.umpari".  n.  u'r.e  ?ft:V.'A  research 
comtustor  is  a rectangular  d..ct  witn  -i  cr.  s.'.  section  - by  3.0  in.  and  at. 
overall  length  of  about  10.0  in.  fl  a.meholder  e:;;p>loyed  for  the  measure- 

ments described  herein  was  a ste.,-.  p.,a;*:r  lO.P'"'  i.n.  thicr.l  containing  an  arra;,- 
of  elgh.t  (8)  hole.s  as  entvr.  in  Fig.  P.  .'r.e  holes  are  0.3'-7  in.  in  diameter 
with  0,75  in.  between  hole  cer.teiu;.  . r.  PbV.'A  o .mtusticn  experiment  employed 
prevaporized  Je„-A  aii'oraft  fuel  which  wau  premixed  with  a preheated  air 
stream.  The  noitiinal  f^eu  comp."! si’ ion  was  Sii.''  weight  percent  carbon  arid  li^.l 
weight  percent  hydrogen  with  a moieculai  weight  ..f  l60.0-,  theref.  re,  the 
nominal  chemical  formula  f i-  the  fuel  is  The  f.^el  heat  of  com.- 

bustion  is  approx imiately  4.3  x 10'  -ule  'kg  (13500  rtu''lbm),  and  a constant 
fuel  specific  heat  (Cp)  ol  2510. -1  Joule  ^kg-'''?;  (0.6  Btu,  Ibn.-'^K)  was  used  in  the 
present  calculations.  The  air  prehea-  temperature  measured  far  upstrean.  of 
the  flame  holder  was  about  756'^K  the  total  .tziss  fl-.w  rate  through 

the  ccm.bustor  was  2.86  x lO"*^  xg/sec;  and  the  fuel  mass  fraction  was  0.0322. 

I'nfortunately , subsequent  tnerm.  c.  uple  measurtsients  ir.  ’ne  PkeVvA  combustor 
nave  showr.  tne  mixture  temperature  direcr.iy  upstrea.m  of  ■’ne  flameholder  tr  be 
approx iTiately  1020°K.  (1400'^F).  The  higher  temperature  of  1030^K  believed  to 
be  present  upstrea.ai  of  the  flameholder  ma-.'  be  attrll  uted  1 hea-*  transfer 
through  the  flameholder  fr-.m  the  com.  us t. ion  z-  rie  1..  the  gas  mixt  ure  upstrca.i. 
of  the  flameholder,  and  p s.sibly  fuel  pyr.  lysis  ups' ream  .f  the  flame- 
holder.  The  drwnetrea.T.  equiiibriun.  temperature  results  tained  previousl'.' 
(Ref.  19)  with  an  inlet  ‘emperature  l'  7^'0^'n  a,-.d  adiatatic  wall  conditior;s 
indicate  that  there  are  signifioatit  er.ero'  losses  in  tl.e  c,  ml  us’i  r.  Propei- 
treatment  of  the  flarr.eholder  bo'a.ndary  c.  ndit  ions  would  require  additional 
information  regarding  the  hea’  transfer  r surface  temporal rre  of  the  flame- 
holder. If  there  are  significant  energy  losses  thr  ug.h  the  flameholder, 
these  would  not  be  correctly  represented  by  tne  udiat  a*.,ic  wall  conditions 
actually  employed  in  the  calculation. . However,  correct  incorporation  of  the 
heat  transfer  mechanism  would  be  sensitive  to  the  local  gar,  temperature,  and 
therefore,  would  require  n reasonable  repi esentat ion  of  the  hydriicarbon 
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comk'US^  ion  Kinetics.  F.>r  complex  IV-lr.  Puch  an  -A  *hir  in  an  extremely 
COTiplica* ‘•'1  problem,  since  the  actual  t'uei  pyrolysis  mechanism  in  generally 
not  known.  Typically,  one  would  be  forced  to  assume  a global  partial  oxida- 
tion react  ion  with  CO  and  H2  as  products.  Then  a detailed  reaction  mechanism 
for  oxidation  of  CC  a.nd  Ho  bo  form  the  fi.nal  products  COg  and  H^O  could  be 
employed.  It  should  be  apparent  tha»  even  sucn  a relatively  "simple"  approach 
to  the  combustion  chemistry  is  still  prciiibitive  in  the  context  of  a elliptic 
fluid  mechanics  calculation  procedure,  altiiO’ogh  some  recent  progress  has  been 
made  in  two-dimensional  modeling  of  a well-stirred  reactor  by  Wormejk  and 
Pratt  (Fef.  1+2''. 


The  computational  region  boundaries  considered  in  the  present  calculations 
are  indicated  by  dashed  lines  in  Fig.  2.  In  order  to  reduce  the  number  of 
grid  points  required  in  the  calculation  the  hole  pattern  was  assumed  t'  be 
periodic  in  the  lor*gitudinal  direction,  so  tha’  syttmetry  b undary  conditions 
could  be  applied  and  only  one -half  of  the  inlet  port  need  by  considered,  as 
showTi  in  Fig,  2,  The  temperature  measurements  in  the  combustor  indicate  that 
this  is  a reasonable  approximation.  The  computational  region  and  coordinate 
system  are  illustrated  in  Fig.  3-  Figure  4 shows  the  sections  (A,  B,  C)  for 
which  profile  plots  of  axial  velocity,  temperature,  and  nitric  oxide  are 
presented.  All  predictions  presented  were  made  using  a 1?  by  9 by  17  grid 
(2601  grid  points^  for  the  x^,  x^,  Xv  directions,  respectively.  A maximum 
axial  velocity  (Uj.)  of  159-5  m/.<:ec  wan  required  to  match  the  exjierimental 
mans  flow  rate.  The  reference  length  (L)  for  all  coordinates  is  0.010935 
meters  (0.75  in.).  The  calculation  wa.n  initiated  downstream  of  the  flame- 
holder  assuming  a mixture  inlet  temperature  of  1025°K  and  a pressure  of 
1.002  X 105  Pascals  (atmospheric  pressure).  A somewhat  arbitral^/  but  plausi- 
ble free  stream  turbulence  kinetic  energy  level  of  ten  percent  (i.e.,  kpg  = 

0.1  ug/2  ) was  imposed  at  the  computational  inlet  plane. 


Axial  velocity  profiles  are  given  in  Figs.  5,  6 and  7 at  sections  A, 

B,  and  C,  respectively,  for  a series  of  axial  stations.  A series  of  cross- 
sectional  plane  contour  plots  of  constant  axial  velocity  are  shewn  in  Figs.  8 
to  11  for  various  axial  stations.  The  qualitative  characteristics  of  the 
Jet  expanding  into  the  combustor  seem  to  be  represented  quite  reasonably. 

As  one  proceeds  away  from  the  inlet,  the  Jet  spreading  is  evident  and  reverse 
flow  is  clearly  present  downstreari  of  the  flameholder  waii  region.  Due  to  the 
large  Inlet  velocity  and  the  relatively  .small  number  of  grid  points  employed 
in  the  calculation,  cell  Reynolds  n’jmbers  In  the  axial  direction  iiecomc  iarge 


L< 


42.  Wormeck,  J.  J.  and  D.  T,  Pratt;  Computer  Modeling  of  Combustion  in  a 
Longwell  Jet-Stirred  Reactor.  Proceedings  Sixteenth  Symposium  (inter- 
national) on  Combustion,  The  Combustion  Institute,  1976,  to  be  published. 


21 


I 


! 

i 

I 


so  that  sirT^i  fl  oat.t  9r;::':.-'ci  p : •r’-rc-r.  <- 

I equations  (see  St"’!_io;_  T1 : • u ■ ' 'ti  •.ne 

present  case  resulted  i:;  a vru'.,  ..I'V,:-  ^ j,*;.  ■ it.  ra,_- 

flux  between  the  t r,let  et..-  c > ■ 'ut  at : (v  ; r , r.cr,  h r.ensiona:  ■ , 

due  to  the  artificial  11  f;' -us  ' on  t'.-r,  rrU  : 'he  rO:r.  .ivy  (:q-uation. 

Profiles  of  the  calculated  turbul-ui.  f i-.:  •■ner..,,  li,  r/r..  iintor  ar*. 

shown  in  Figs.  12,  13  and  . Th*'  ' :o'l  ut.vicc-  of  t:...  a,  , . In'.  <;t  free 

strean:  turbulence  level  (!?;•■'  re<on.-  n yernlnt  fa’.  :■  • ;;v,7;,cnrcar.  of  tru- 

inlet.  The  effect  of  kgg  vaj'iatlnri.  ..  ‘,h»-  pertln>ei;  r.','.  , . en  variables 

was  not  investigated  .In.  the  p.p-f.'  on'  ;n.  .dy,  :;u‘.  tiil..  q --;-Tit,r.  clearly  neea.c 
to  be  addrttssed  in  fi'.ture  experinrx-nta . a:.u  t.-.eoretxca..  studies . 

The  equilibrium  hydrocarbon  cheraistry  asEuinxjtion  van  not  expected  to  be 
valid  near  the  inlet  port  of  the  cotnbustor  under  consideration.  Therefore, 
an  ad  hoc  ignition  delay  criterion  v/as  imposed  or.  the  solution  to  prevent 
unrealistically  large  temperature  gi-ariients  from  occurring  near  the  inlet. 

This  was  accomplished  in  the  present  numerical  procedure  by  specifying  the 
pseudo-kinetic  rate  constant,  (Rei'.  I9),  as  a function  of  axial  distance 
from  the  inlet  plane.  The  effect  of  *he  specified  ignition  delay  on  the  port- 
centerline  axial  temperature  vari atio-ns  is  showr.  In  Fig.  It,  along  with  exper- 
imental measurements  at  four  axia.1  stanlons.  As  e.xpected  the  predicted  tem- 
peratures are  considerably  higher  than  the  meaHured  values  because  ol'  tt.e 
energy  losses  in  the  combustor  whic’n  were  not  Included  in  the  calculations. 
Nondlmenslonal  temperature  profiles  are  presented  in  Figs.  I6  to  l6  at 
actions  A,  B,  and  C,  respectively,  for  a series  of  axial  stations,  and  a 
sequence  of  cross-sectional  plane  contour  plots  of  constant  temperature  ore 
shown  in  Figs.  I9  to  21.  Tiie  availebj.e  experimenta.i,  teniperature  measurements 
are  also  displayed  on  Fig.  17  for  section  B.  The  presence  of  artificial  vis- 
cosity (in  the  axial  direction)  in  both  the  .species  and  energ;/  equations  may 
have  distorted  the  transverse  tem.perature  profiles,  since  outside  the  iet 
boundaries  artificial  diffusion  is  siiiyiif icar.t  compared  to  convection.  How- 
ever, the  qualitative  comparisons  which  con  be  mode  frar.  r'ig.  1?  Indicate  a 
reasonable  variation  of  predicted  temperature  across  the  combustor.  Only  a 
limited  niunber  of  nitric  oxide  (IJO)  concentratl  on  mea-suremonts  were  made 
in  the  combustor  under  conrui deration.  Since  the  NO  co.uientration  is  sensitive 
the  temperature  history  of  a fluid  particl»>,  comparl .sons  between  predic- 
tion and  exjjeriment  are  not  .meaningful  in  this  in.stauce.  The  NO  concentration 
profiles  (Fig.  22)  are  sho-wn  only  for  tlio  axial  station  where  an  equilibr.imn 
temperature  has  been  aciiieveci. 
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A. 


The  present  numerical  results  aemonstrat.e  both  the  basic  integrity  oi'  the 
computat lonal  procedure  developed  under  the  pr-esent  study  and  ’he  capability 
of  the  MINT  code  to  perform  calculations  of  three-dimiens i onal  combusting  flows 
with  recirculation.  The  need  foi-  improvement  of  the  differencing  scheme  is 
apparent  because  of  the  suspected  distortion  in  the  predictions  due  to  exces- 
sive artificial  viscosity.  The  preliminary  efforts  to  develop  a simple  com- 
pact differencing  scheme  made  under  the  present  progra.m  showed  some  success  in 
nonreacting  flows,  however,  no  reasonable  results  could  be  obtained  for  com- 
busting flows.  The  ad  hoc  nature  of  the  compact  difference  scheme  proposed  in 
Section  III  is  an  undesirable  feature,  and  it  may  preclude  successful  use  of 
schemes  such  as  these  for  general  problems.  The  motivation  for  continued 
development  of  higher  order  solution  procedures  should  still  be  high  for 
the  reasons  mentioned  in  Section  III. 

In  addition,  there  is  clearly  a need  for  an  improved  hydrocarbori  chemistry 
model  in  the  existing  computer  code  if  reasonably  accurate  quantitative  predic- 
tions are  desired.  For  a general  three-dimensional  elliptic  computational 
procedure  such  as  the  MINT  code,  it  is  felt  that  a relatively  simple  finite- 
rate  hydrocarbon  chemistry  model  should  be  employed,  in  conjungtion  with  a 
simple  model  for  representation  of  the  influence  of  torb’ulent  concentra’ ion 
fluctuations,  such  as  that  employed  by  Hutchinson,  JOialil  and  Whitelaw  (Ref.  U3) 


U3.  Hutchinson,  P.,  E.  E.  Khalil  and  J.  H.  Whitelavr;  The  Calculation  of  Wall 
Heat  Transfer  Rate  and  Pollution  Formatioti  in  A.xisymmetric  Furnaces. 
Presented  at  4th  Members  Conference,  International  Flame  Research 
Foundation,  Ijmuiden,  Holland,  Kay  1976. 
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FREP  COMF,rrER  CGDF 


The  two-dimensional  FREF  (Field  Relaxation  Klliptic  i'rocedure ) 
computer  code  had  been  developed  at  United  Technolofrie.;  Research  Center 
under  both  Air  Force  and  Enviroamental  In:otection  Af^'enoy  sponsorship  for 
the  prediction  of  combusting  flow  fields.  Presently,  the  FHEP  code 
is  being  further  developed  under  EPA  contract  nwbor  68-02-1873  and  this 
code  is  deliverable  to  the  Air  Force  'under  tne  present  contract.  The 
theoretical  analysis  for  the  FREP  code,  comparison  of  comp\ited  results 
with  measurements,  and  tne  user's  manual  for  the  FI{EP  comjiuter  program 
will  be  available  as  EPA  reports. 

In  order  to  demonstrate  t.he  fRF.P  comi-uter  co>le  on  the  Wright- 
Patterson  Air  Force  Base  CDC  6fi00  a sample  two-dimensional  problem  was 
rmin.  TTie  ge'-metry  considered  was  similar  to  that  in  tiie  tnree-dimcnsional 
computation  with  each  array  oi'  four  holes  replaced  by  a silt,  and  the  end- 
walls  -were  neglected.  The  inlet  temperature  was  taken  as  7%°K  (900*^F)  witn 
all  other  conditions  as  described  in  Section  IV.  An  arbitrary  linear 
ignition  delay  criterion  (from  inlet  to  exit  plane)  was  applied  in  con, junc- 
tion with  equilibrium  hyurocarbon  cliemiutry  for  this  deiaonstration  calcu- 
lation. The  computed  axial  velocity  profile.")  at  several  axial  stations 
are  shown  in  Fig.  23;  the  plane  of  symmetry  and  the  channel  wall  are 
located  at  Xj  - 0.0  ana  X_,  1.0,  respectively,  while  tne  inlet  port  is 

centered  between  these  two  boundaries.  Tne  qualitative  behavior  of  these 
profiles  seems  reasonable,  in  view  of  tne  coarse  cfxnputational  mesh 
(17x17)  employed  in  this  calculation.  Extensive  evaluation  of  the  FREP 
code  is  being  carried  out  under  EPA  contract  68-02-1873. 
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Figure  2.  Experimental  ftameholder  configuration  for  three  dimensional 
rectangular  combustor 
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SECTION  C DIAGONAL 


(Mondimensional  axial  velocity  profiles 
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Figure  12.  Nondimensional  turbulence  kinetic  energy  profiles 
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Figure  14.  Nondimensional  turbulence  kinetic  energy  profiles 
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Figure  17.  Nondimensional  temperature  profiles 
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Figure  21.  Contour*  of 


I 


SECTION  B X,  - 0,5 
— X3  4 0,  PREDICTION 

0 I 

□ X3  - 4,0  ' EXPERIMENT 

A X3*  5 33  I 


Figure  22.  Comparison  between  predicted  and  experimental 
nitric  oxide  (NO)  concentration  profiles 
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Figure  23.  Nondimensional  Axial  Velocity  Profiles  for  Two-Dimensional  FREP  Calculation 
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